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Abstract We present a new method to automatically track filaments over the 
solar disk. The filaments are first detected on Meudon Spectroheliograph Ha 
images of the Sun, applying the technique developed by |Fuller, Aboudarham, and] 



Bentley {Solar phys. 227, 61, 2005). This technique combines cleaning processes, 
image segmentation based on region growing, and morphological parameter ex- 
traction, including the determination of filament skeletons. The coordinates of 
the skeleton pixels, given in a heliocentric system, are then converted to a more 
appropriate reference frame that follows the rotation of the Sun surface. In such a 
frame, a co-rotating filament is always located around the same position, and its 
skeletons (extracted from each image) are thus spatially close, forming a group of 
adjacent features. In a third step, the shape of each skeleton is compared with its 
neighbours using a curve-matching algorithm. This step will permit us to define 
the probability [P] that two close filaments in the co-rotating frame are actually 
the same one observed on two different images. At the end, the pairs of features, 
for which the corresponding probability is greater than a threshold value, are 
associated using tracking identification indexes. On a representative sample of 
filaments, the good agreement between automated and manual tracking confirms 
the reliability of the technique to be applied on large data sets. Especially, this 
code is already used in the framework of the Heliophysics Integrated Observatory 
(HELIO) to populate a catalogue dedicated to solar and heliosphcric features 
(HFC). An extension of this method to others filament observations, and possibly 
the sunspots, faculae, and coronal holes tracking can be also envisaged. 

Keywords: Solar filaments, Halpha observations, automated tracking, image 
processing, virtual observatory, HELIO, HFC 
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1. Introduction 



The heliophysics community embraces a number of existing disciplines - so- 
lar and heliospheric physics, magnetospheric and ionospheric physics for the 
Earth and other planets - and thus must deal with very large sets of het- 
erogeneous data from ground- and space-based instruments. The Heliophysics 
Integrated Observatory (HELIO: http: / /www.helio-vo.eu) is a virtual observatory 

HELIO 



2011) 



dedicated to the solar physics and heliophysics ( Bentley et al. 
provides a distributed network of services, which helps the researchers to eas- 
ily mine relevant information and data. Especially, the Heliospheric Feature 
Catalogue (HFC: [http:// voparis- hello. obspm.fr/hfc-gui/index.php ) is a database- 
oriented service that allows access to a large amount of solar and heliospheric 
features data. Extraction of feature information stored in the HFC is realized 



using an increasing number of recognition codes ( Fuller, Aboudarham, and Bent 



ley 




2005 


Zharkov et al. 2005 Barra et al. \ |2009 Krista and Gallagher 


Lobzin et al. 




2009 Higgins et al. 




20111. 



We present here a new algorithm for the solar-filament tracking. It has been 
initially developed in the framework of the HELIO project in order to provide 
tracking data, as a supplementary information, to the description of filaments 
already available in the HFC. The filaments are large-scale structures of relatively 
dense and cool plasma suspended in the hot and thin corona. They are particu- 
larly visible on Ha observations, where they appear as elongated dark features 



with several barbs on the solar chromosphere ( Tandberg-Hanssen[ 1995 Mackay 



et al. 20101. During their lifetime, their shape and intensity can be subject to 
a number of modifications, especially, part or all of a filament can sometimes 
undergo sudden disappearances (which may be followed in some cases by re- 
appearances) at some wavelengths, probably due to changes in their background 
characteristics - temperature, pressure, etc. - or caused by eruptive process. 
Hence, this sudden and unpredictable behaviour can make filaments difficult 
to track over the Sun surface. However, following these features over time is 
relevant, notably for space weather, since erupting disappearances can play a 



role in the triggering of coronal mass ejections (CMEs) (Gilbert et al. 
Gopalswamy et al. 20031. 



2000 



Many automated methods have been successfully developed during the last 



ten years to detect filaments, including Gao, Wang, and Zhou 



Aboudarham, and Bentley ( 2005 ), Bernasconi, Rust, and Hakim '( 2005|),|Zharkovc | 



(2002L Fuller, 



Rust, and Hakim (2005); Joshi, Srivastava, and Mathew (2010) also propose 



and Schetinin (2005), and more recently Joshi, Srivastava, and Mathew (2010). 
In addition to their recognition algorithms, Gao, Wang, and Zhou (|2002j ); |Bernascj oni 



tracking capabilities that allow to identify disappearances by following feature 



locations day by day. In particular, Gao, Wang, and Zhou (2002) have tracked 
filament disappearances over a full year, and the code of Bernasconi, Rust] 



and Hakim (2005) is currently applied to detect and track filaments on Big 
Bear Solar Observatory (BBSO) images, providing among others data to the 
Heliophysics Events Knowledgebase (HEK: http://lmsal.com/hek/). In all cases 
the basics of the methods are quite similar: starting from heliocentric positions of 
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the filaments detected at a given time (i.e., on a given image), they estimate the 
coordinates of their respective centroids on the previous and/or next adjacent 
images using equations that correct for the solar rotation ( Cox[ 2000). Then, 
they check if detected filament lies within a given circle area from the predicted 
locations; if it does, the two filaments are considered to be the same. In the 



case of Bernasconi, Rust, and Hakim (20051, where no filament is found they 



then extend the search up to three days to confirm that the filament actually 
disappears or not. As highlighted by the authors, the 3-day search is motivated 
by the fact that filaments sometimes change shape so much between times of 
observation, that their location may fall out of the search area, temporarily 
loosing the tracking. A consequence is that during this period of time it is not 
possible to know if the filament actually disappears, or if its shape is just strongly 
deformed due to splitting or partial disappearance. A tracking algorithm based 
on time, position, but also shape matching would permit to distinguish between 
the different possible behaviour. This capability appears to be essential to detect, 
among others, erupting disappearance. 



In our method, filaments are first detected on Ha images by an automated 
recognition code that extracts, among other things, the pruned skeletons of fila- 
ments. A pruned skeleton can be defined as the thickest curve inside the feature 
area that preserves the shape topology; in terms of image processing it results of 



a thinning followed by a pruning operation ( Gonzalez and Woods, ,2002) applied 
on the feature pixels. Then, the coordinates of the resulting skeleton pixels are 
converted to a more appropriate reference frame that rotates with the solar 
surface. In such a frame, the skeleton of a filament observed on successive images 



appears as a cluster of close curves (Mouradian 1998). Hence, a comparison 



in this frame of each skeleton curve with its closer neighbours using a curve- 
matching algorithm, allows us to calculate probabilities that all of these features 
actually belong to the filament. The algorithm demonstrates its ability to track 
filaments on large data sets, limiting the number of false tracking detections. 
Section [2] will introduce the observations used to test the technique, and briefly 
describe the filament detection process; Section [3] is devoted to the explanation 
of the tracking algorithm; Section [4] will present the resulting performances of 
the code. Finally applications will be discussed in Section [5] 



2. Detection of Filaments 
2.1. Observations 

The Meudon spectroheliograph of the Observatoire de Paris performs daily ob- 
servations of the solar photosphere and chromosphere at three wavelengths: blue 
wing of the Can K 393.4 nm line, or Ki„ ; centre of the Can K 393.4 nm line, 
or K3 ; and Ha 656.3 nm. More than fifteen years of data are accessible from 
the BASS2000 website (http: / /bass2000. obspm.fr/home.php), which permits us 
to apply the tracking algorithm over a large period of time. 
The images produced are 2D heliocentric projections of the full solar disk as 
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Figure 1. a. Gray-scale image of the Sun observed in Ha by the Meudon spectroheliograph 
on 21 August 2002 at 11:59:00 UT. The image is centred on the solar disk, and modified to 
have X- and j/-axis respectively aligned with the Equator and rotation axis. Due to their cooler 
temperature, the filaments are seen like dark features on the solar surface. Centre-of-limb 
effects as well as dark faint lines, caused by dust particles on the optics, are also visible, b. 
Same image after cleaning processes. 



seen from the Earth, with a typical size of 1024 by 1024 pixels, and a spatial 
resolution of w 2.28 arcsecs. The origin of the reference frame corresponds to 
the disk centre, x-axis is aligned with the Sun's Equator and points towards 
the west limb, and y-axis is aligned with the rotation axis and points towards 
the North. Figure [I] shows such an image before (left panel) and after (right 
panel) cleaning processes ( Fuller, Aboudarham, and Bentley[ 2005) ; these pre- 
processes are required to avoid false detections, but also to optimize the efficiency 
of the region-growing process used by the detection algorithm, by increasing the 
contrast between the filaments and quiet-sun intensities. 



2.2. Recognition Code 



In order to retrieve the location and morphology of filament skeletons required 
to perform tracking, an automated recognition code is first run over the full 
data set available. We use here an algorithm developed in the framework of the 
European Grid of Solar Observations (EGSO: http://www.egso.org) project, and 
successfully applied on Ha data from the Meudon Spectroheliograph and BBSO 
( Fuller, Aboudarham, and Bentley] 2005). This code is now also used in HELIO 
to provide to the HFC, a description of the filaments detected on Meudon images. 
(Filaments detected on BBSO observations are planned to be added in HFC in 
few months.) 

The recognition method requires three main steps to be completed: 

i) Since images are sometimes blurred, which can reduce the efficiency of the 



filament segmentation, a Laplacian spatial filter ( Russ 2002 1 is first used to 
enhance the clearness of filament contours. 



ii) Then, a region-growing algorithm (Gonzalez and Woods 2002) permits us to 



group pixels of a same filament together. Starting from a seed region (found 
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using an appropriate threshold value), the procedure searches the connected 
pixels for which the intensities are within a range defined by the mean and 
standard deviation of their neighbours. 



Hi) A morphological-closing operator (Gonzalez and Woods 2002) will finally 



merge resulting nearby regions that could be considered as a single filament 
on the segmented images. 

From the boundary detection, several parameters can be extracted concerning 
the intensity, location, or morphology of filaments. In particular, the algorithm 
uses thinning and pruning methods (Gonzalez and Woods 2002) in order to 



compute skeleton shapes (results of a skeleton computation is illustrated on 
Figure [2|. Two conditions are required to produce such shapes: all skeleton 
parts should be connected, and the full extent of the region must be contained 
in the skeletonized representation. Knowing this characteristic will permit us to 
define, among others, the length, the centre, and the curvature of the filaments. 
At the end, all of these parameters will be ingested in the HFC from where they 
can be downloaded using dedicated query interfaces. 

For our purpose, the shape of the skeletons will serve as a matching criterion to 
identify a filament from an image to the following one. 



3. Tracking of Filaments 
3.1. Description of the Algorithm 

In this section we present in detail the tracking algorithm. In addition to the 
coordinates of skeleton pixels provided by the recognition code, two types of 
indices are also required to follow co-rotating features over time. The first index, 
called feature identification number [I'i], is allocated to each filament [i] extracted 
by the recognition code. All of its values are unique (the values of the feature 
identification number follow typically the order of detection of the filaments) in 
order to identify and retrieve in the HFC, a feature detected at a given time 
on a given image. The second index, called tracking identification number [r,], 
is assigned to each filament [i] by the tracking code. As opposed to i^i several 
filaments can have the same value for t^, which will be an indication that these 
filaments are actually components of the same co-rotating feature observed on 
successive images. The values of are set by the tracking code, comparing and 



associating filaments two by two (as explained in Section 3.1.4). At the beginning 
of the execution, all of the tracking indexes [r^] are initially equal to the feature 
ones [i'i] by default. 

The main steps of the tracking method can be summarized as follows: 

i) We select the time [t^] of the image [l] for which we want to perform tracking. 

a) We then define a time window [ioj^i] centred on and spanning one Car- 
rington rotation of the Sun, such that to = ti. ~ 0.5Tcarr and ti=tc + 0.5Tcarr, 
where Tcarr = 27.2753 days. The choice of this time range makes sure that all 
of the filaments detected at will be fully tracked over a period of O.ST^arr 
days (which is actually the average time required to cross the solar disk). 
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Figure 2. a. Same image as in Figure but over-plotted with coloured contours and skele- 
tons extracted from the recognition code. (The colors used here are randomly distributed to 
distinguish individual features detected.) b. Zoom on the area marked by the black box on the 
image displayed in Panel a. 
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Hi) We download the HFC information about filaments detected between to and 
ti] parameters required are typically the pixel coordinates of the pruned 
skeletons, but also their corresponding feature identification numbers [v]. (In 
practice, the corresponding tracking identification numbers [r] are also loaded 
to ensure a continuity with possible previous tracking runs.) 

iv) Since the pixel coordinates are given in the reference frame of images (see 
Section 2.1 1, we convert them to a more appropriate coordinates system as 
explained in the following section. 

v) We apply a curve-matching algorithm in the new frame. The purpose of this 
algorithm is to compare the shape of all of the possible pairs of skeletons: if 
two skeletons have similar shapes, then they are assumed to be two distinct 
parts of the same filament moving over the Sun's surface. 

vi) When the comparison succeeds, we associate the two corresponding skeletons 
by allocating them the same tracking identification number. 

vii) The same process is then repeated on the next image [t + 1] taken at time 
saving all of the tracking numbers [r] defined during previous runs on 
the image [l]. This operation allows us to conserve tracking information from 
a processed image to the following. 



3.1.1. Reference Frame used for Tracking 



The tracking for filaments detected at time t — is performed on a specific 
reference frame 3? that follows the solar rotation between and ti . Such a frame 
has the advantage to: i) stabilize the coordinates of co-rotating filaments around 
the same position, by correcting the translation of heliocentric longitudes for the 
rotation, ii) make appear all the segments of a filament visible between tp a-nd 
ti, which may suffer full or partial disappearance during its lifetime, iii) since 
along the first dimension time and space are intermingled, offer fast computation 
time by working in two dimensions instead of three. (The computation for all of 
the filaments detected over one solar rotation takes less than w 30 seconds on 
average on a 3.06 GHz Intel duo core machine.) 

For each image processed the detection code returns the heliocentric coordinates 



{X,Y) of skeleton pixels (as defined in Section 2.1 1, which must be converted 
in the proper coordinates system associated to 3? before being usable by the 
matching algorithm. To achieve this goal, the heliographic longitudes and lati- 
tudes {iph^^h) of skeleton pixels on the Sun's surface are first computed using 
an appropriate transformation matrix ( Thompson 2006 ) , then the coordinates 
((p, A) in the co-rotating frame 3? are deduced using the relations: 



A = A^, (1) 

= ^^carr(io " *) + (V'to " "^h) + VORci^h), (2) 

where l^carr is the Carrington rotation speed in degrees.day~^ , t is the time 
when the filament was detected in decimal days (i.e., the time of the observation), 
to is the start time of the tracking period in decimal days, (pto = 360° is the longi- 
tude of the central meridian in 3? at ioi and fDRci^h) is a term that corrects the 
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effects caused by differential rotation (Ulrich and Boyden 20061. These effects 
are not significant here, but will become important when the algorithm will 
be run to track filaments from one solar rotation to the following (see Section 
3.2). We note that the coordinates system defined here is quite similar to the 



Carrington one, but starting at (^t^. 



The steps of calculation of the skeleton coordinates in the reference frame 3? 
are illustrated on Figure [3j In this example, we want to track filaments detected 
on a Meudon spectroheliogram on 21 August 2002 at 11:59:00 UT; date that 
defines here the time of observation [t^]. Thus, according to the discussion in 



Section 3.1 the tracking process concerns all of the filaments extracted between 



7 August 2002 at 20:40:47 UT and 4 September 2002 at 03:17:12 UT, which 
correspond to the start time [to] and the end time [ti] of the tracking period 
respectively. Panels a, b, c, d, and e, show a sample of five images of the Sun 
obtained during this time range; the date of observation is indicated in the 
upper-left corner of each image (for more clarity only a few images produced 
during this period are displayed). The heliocentric coordinates of skeleton pixels 
as well as the identification numbers of each detected filament are also plotted 
using a coloured scale. Panel f presents an example of corresponding coordinates 
((p, A), calculated for the co-rotating filament visible inside the white bounding 
rectangles on the previous panels. In 3?, all of the skeletons belonging to this 
filament form a structured cluster around the same location, giving the full 
shape of the filament over the disk crossing. 

Figure |4] presents the resulting synoptic map of skeleton pixel's coordinates 
converted in 3? for the previous example over the full period [io,ii]. The lon- 
gitude axis uses the Carrington convention starting at Lptg = 360°, and ending 
at ifti = 0°, with (ft^ = 180°, corresponding to the longitude of the central 
meridian on 21 August 2002 at 11:59:00 UT. At this stage of the procedure, 
the heliocentric and co-rotating pixel coordinates of features, as well as their 
corresponding identification numbers [v and r], are stored by the tracking code. 



3.1.2. Curvilinear Interpolation 



To be consistent, the curve-matching algorithm must be applied on two curve 
segments having the same number of points and the same length. However, due at 
least to the coordinates' projections, the distribution of points along the skeleton 
curves is not necessarily uniform, hence, we need to interpolate the curves in such 
a way that the distance [ds\ between two consecutive points is constant. At the 
same time, such a procedure can be seen as a smoothing operation that will 
highlight the main shape of the skeletons. It results that the value in degrees of 
[ds] to use in the co-rotating frame has to be a compromise between the spatial 
resolution (i.e., minimum distance between two points) and the length of the 
skeletons: too small a value can significantly increase the computation time, but 
also the possible spatial fluctuations along the curve at small scales, whereas too 
large a value decreases the number of points along the curve, and so the spatial 
resolution used to define its shape. 

Once an acceptable value of [ds] is found, starting from the location of the first 
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|2002-08-25T07:23:00 




Figure 3. Illustration of the calculation of the skeleton pixel's coordinates from the heliocen- 
tric coordinates system to the co-rotating one. Panels a, b, c, d, and e show a sample of five 
images of the Sun produced between 7 August 2002 at 20:40:47 UT and 4 September 2002 
at 03:17:12 UT (on each image the date of observation is written in the upper-left corner). 
Example of skeleton pixel's coordinates [ip, A) calculated in the co-rotating frame 5R are plotted 
on Panel f; the region displayed corresponds to the white bounding rectangles on each image. 
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Figure 4. Synoptic map of tlie slceleton pixel's coordinates in 5R between 07 August 2002 at 
20:40:47 UT (ie., to) and 04 September 2002 at 03:17:12 UT (i.e., ti). Carrington convention is 
used to define the longitude axis. The vertical darker line at the centre indicates the longitude 
•ft^ = 180°, which corresponds the longitude of the central meridian on 21 August 2002 at 
11:59:00 UT. The bounding rectangle delimited by black dashed line shows the area displayed 
on the panel f in Figure |3] 



point Ai) at one of the ends of the curve (the order of the points along the 
curve goes from one end to the other), the interpolation method consists in the 
determination of the intersection point {fki„tJ ^ki^t) between the circle of radius 
[ds] centred on (ipi^Xi), and the line passing through the two points {ipk,Xk) 
and {ipk+i, Afe+i) of the curve, which satisfy the condition: 

k fe+1 

J2 v/(Vi-<^i-i)^ + (A;-A(_i)2 < ds < ^ v/(¥'/-W-i)^ + (Ai-A,_i)2. 
1=2 1=2 

(3) 

When the intersection point (<y5fei„t7 ^ki^t) is known, the process is then repeated 
replacing ((/ji,Ai) by ((ySfcint j ''^fci„t) the new circle centre, and calculating the 
next intersection point {'^ki^i+i ' -^fci„t+i) with the curve. The interpolation finally 
stops when the other end of the curve is reached. 

At the same time, the lengths of the skeletons are estimated by simply summing 
the distances [ds] found between each interpolated points (i.e., ((Pfcint j -^fcint) ^'^^ 
(</5fei„t+n '^fcint+i))- This technique will slightly under-estimated the actual length, 
which appears to be not really significant on the tracking results. 
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3.1.3. Curve Matching Algorithm 

In order to significantly reduce the computation time, the matching algorithm 
is only applied on skeletons that are close enough on the reference frame. To 
achieve this goal, for each coordinates ((^fc^, Afc. ) of a first skeleton i, we search 
the coordinates {ipk, , ^k, ) of others skeletons j {j 7^ i) for which the distance 
Tk.kj = \/iVk, - VkjY + (Afc, - XkjY is less than a maximum value r^ax (de- 
fined in degrees in 3?). If the condition r^.^^ < r^ax is fulfilled, then the curve 
matching algorithm is applied. In practice the value of the input parameter 
''max does not affect significantly the efficiency of the algorithm, because it will 
automatically dissociate features that are too much distant. 

Given two curves Ci and C2 of close skeletons 1 and 2 from the same fila- 
ment observed on two different images at two different times, there should be a 
Euclidean transformation E, such that ECi matches C2. However, since during 
its lifetime a filament can be subject to shape deformation, segmentation, but 
also translation and/or rotation, no exact match is likely to occur, and so we 
look for a transformation E giving the best match in the least-squares sense. 
Let Ci and C2 be represented by the vectorial sequences (uk)^^i and (vk)^^i 
respectively, where n is the number of points of the curves (we assume here that 
the two curves have the same length) . Matching consists of finding a Euclidean 
transformation E of the plane that will minimize the distance [P] between the 
sequences {E\\\^)^^^ and (vk)fc^i, such that: 



A = min^Z^ = min^ ^ jE'uk - Vkp. (4) 

fc=i 

To simplify the calculation, the curve Ci is first translated such as X]fc=i — 0. 
Then, if we write E as i?Uk = i?eUk -f a, where Rg denotes a counter-clockwise 
rotation by an angle 9, and a a translation; the best match is obtained when 



Ayache and Faugeras 


(1986 


); 


Wolfson 


(1990 



a= - V'vk, and 6* = - V'ufeU/s, (5) 
n ^-^ ^-^ 
fc=i fe=i 

where Ufe and Uk are complex representations of the vectors Uk and Vk respec- 
tively. 

At the end, for each pair (C.j, Cj) of skeleton curves analysed, the matching 
process returns three best-fitted parameters 

(a=|a|+(5a,0 = |0M = \/A) , (6) 

where (5a is a term that corrects the 'Y^!k=\ = translation introduced previ- 
ously. These parameters, which represent respectively the degrees of translation, 
rotation, and deformation between the two curves, will be used to compute in a 
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second step the probability [Pij] that the skeletons [i and j] belong to the same 
filament. 

If two skeletons have different lengths, k ^ Ij, the matching is also performed 
but between the smallest of the two, say i, and all of the possible sub-segments 
k of the other one j that have the same length 1^ = U- In this case, only the 
probability [Pik] for the sub-segment returning the best match is finally retained. 
This technique becomes actually consistent only if the lengths of the two curves 
do not differ significantly, we therefore introduce an additional condition to the 
minimum on maximum lengths ratio: if this ratio is too small, then the matching 
algorithm is not applied. 



3.1.4- Confidence of Tracking 

For each set of best-fitted parameters (a, 9, we calculate a normalized prob- 
ability [Pij] defined by: 

Pzj^WaPa + WePe + WiPi (0<P,,<1), (7) 

where Pa, Pe, and Pi are the probability terms respectively related to the trans- 
lation [a], rotation [9], and deformation [I] parameters, and Wa, We, and Wi 
are the corresponding probability weights, which satisfy Wa + Wg + Wi = 1. 
Pij can be seen as the level of trust of the tracking; the larger it is, the more 
the automated code is confident about the fact the two skeletons [i and j] are 
actually two segments of the same filament. 

Since no formal relation exists, the Pa, Pe, and Pi terms are simply described here 
using linear expressions, nevertheless, under the condition that the probability 
[Pij] decreases when the three parameters a, 9, and I increase: 

= 1 - — (I - Pe) if a < °° , and P^ = otherwise, (8) 

Q 

Pe = 1 - — (1 - Pe) if 6* < , and Pe = otherwise, (9) 



Co 


1 


-Pb' 




Oo 


1 


-Pe 




k 


1 


~Pb 



70 

P, = 1 - -^(I - Pe) if Z < , and P; = otherwise, (10) 

'0 

where oq, 9o, and Iq are input parameters that satisfy Pa=ao = Pe, Pe=ea ~ Pe, 
and Pi^ig = Pe) respectively. At the same time, the three conditions a < jz^, 
6 < i^p^ , and I < will ensure that Pa, Pe, and P; have always positive 

values. The determination of these three parameters' values will be explained 
later in Section [472) Concerning Pq, we take Pe — 0.5 which is actually also the 
value of the respective probability [Pij] above which two skeletons are assumed 
to match, and are thus associated. 

The association step is very important since it will permit us to follow a filament 
on successive images, assigning the same tracking identification number to its 
skeletons. An example of such a process is described on Figure [5j Where two 
skeletons [i and j] are not associated together {i.e., ti ^ Tj), but both are 
associated with a third one [k] [i.e., Ti — and Tj = r^), then all three will be 
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Success 

(Pii i 0-5) 



Ti = min(Ti, Tj) 
Tj = min(Ti, Tj) 





Ti = i 




Tj = j 



Figure 5. Given two skeletons \i and j] of feature and tracking identification numbers initially 
equal to {vi = i,Ti = i) and (uj = j, Tj = j) respectively, where i ^ j, the association consists 
of assigning the same tracking identification number to both skeletons when the matching 
succeeds. In practice the minimum value between Ti and tj is conserved. In the case where at 
least one of the skeletons compared has already received a different tracking number, from a 
previous matching with another skeleton k for instance, Ti = t^ ^ i, then we assign to all the 
features [i, j, and k] the same number Ti = Tj = tj. = min(Ti, Tj , r(j). 



related, and will obtain the same tracking identification number (i.e., = Tj = 
Tk). This condition appears to be very efficient for merging several segments of 
a filament that suffers significant shape modifications. 

Once all of the associations between the pairs of curves are done, for each 
skeleton [i] we finally calculate the average value 



<P^ > = 




where Nj is the number of matchings realized between the skeleton [i] and the 
skeletons [j] {j =/= i) for which tj — Ti. This value will be used as an average 
tracking confidence level for the corresponding filament [i]. 

3.2. Tracking Filaments over Two Successive Solar Rotations 

The algorithm can also be adapted to track filaments that perform two consec- 
utive crossings on the solar disk. If we assume that a given filament is observed 
around a time it, but also around — Team then we can extend the previous 
time range [toi^i] such that w — LST^arr and ti w + CST^arr, where the 
average coordinates of the filament skeletons in the current and in the previous 
rotations will be respectively (<(/?>, < A>) and {<ip> +360°, < A>) (using the 
reference convention defined by Equation ([I]) and Equation ([2])). Figure |6] shows 
the skeleton pixel's coordinates in 5R between to and ti, taken at 11:59:00 UT 
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2002-09-04103:17:12 t t-Tcarr 2002-07-11114:04:21 
— \ — ' ' 1 — I — P — \ — ' ' 1 P ' ' — I — I ' ' •- 




Longitude (degrees) 

Figure 6. Synoptic map of the new skeleton pixel's coordinates between tg at 14:04:21 UT 
on 11 July 2002, and ti at 03:17:12 UT on 4 September 2002. Carrington convention is still 
used here to define the longitude axis. The vertical red line on the left indicates the longitude 
(fit = —180°, which corresponds the longitude of the central meridian at t, on 21 August 2002 
at 11:59:00 UT. The vertical blue line on the right indicates the longitude Vt— Tcarr ~ -|-180°, 
which coincides with the longitude of the central meridian at i — Tcarr- The two groups of 
skeletons coloured respectively in purple and orange belong to the same filament observed 
over two solar rotations (all the others skeletons are plotted here in green for more clarity). 
In this example, the matching process is first used to group together the purple skeletons on 
one side, and the orange skeletons on the other. Then, when each group has its own tracking 
identification nu mb er, the matching algorithm is run to compare the curves of the two clusters 
using Equation | |14| l. In this case the matching process can be done because both group gravity 
centres are separated by 360° . 

on 21 August 2002. 

At this stage, we first apply the matching algorithm in order to set the tracking 
identification numbers (as described on previous sections) , then from a first group 
[gi] of skeletons [i] having the same tracking identification number — Tg. and 
located around {<ipi> , < Xi>), we look for a possible group [gj] of skeletons [j], 
satisfying tj = Tg^ and situated around the coordinates (< (pj >«< (pi > +360°, < 
Xj >«< Ai >). If such a group [gj] is found, then the matching algorithm is run 
between the skeletons of both groups only. In practice, to start the matching 
process, the average gravity center (< ipj >, < >) of the group [gj] must be 
inside the circle of radius rmax and centred on {<ipi> +360°, < A^ >). 

However, since the chance to observe the same filament over several successive 
rotation periods decreases with time, we introduce here a new probability term 
[-Pa*] that is proportional to the difference At of the observation times between 
the two skeletons to match. (Since one solar rotation period separates the two 
features. At should be normally around Tcarr-) Hence, the resulting probability 
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[Plj] that we calculate for each pair of skeletons of the two groups will be 
slightly different: 

Plj = PAtP^J, (12) 

where Pij is the probability defined in Equation ([?]). 

In order to reproduce the filament's average lifetime, the term Pa* must be 
maximal (i.e., Pa* = 1) when At = 0, and must tend towards zero when At 
increases. As for Pa, Pg, and Pi (see Equation (|8])), we use here a linear expression 
defined by: 

Pa* = 1 - ^(1 - Pe) if At < , and Pa* = otherwise, (13) 

Ato 1 - Pe 

where A^o is an input parameter above which the probability [Pa*] falls below a 
given value Pq = 0.5. The value of Ato is also estimated empirically as explained 
in the next section. 

Finally, the code will associate the two groups [gi and gj] of skeletons only if the 
average value 

<p'>= — yyp',, (14) 

* 3 i=l j = l 

where Ni and Nj are respectively the number of skeletons in the first and second 
group, is greater then 0.5. In this case, all of the skeletons [i] of the first group 
[gi] will receive a third identification number [uji], which will be equal to the 
tracking identification number of the group [gj] such as u)i = Tj — Tg.. 
At the end of the tracking computation, a filament will be specified by three 
identification numbers [vi, Ti, and uji] to be identified on a image, to be tracked 
over the solar disk, and to be possibly linked to another feature on the previous 
rotation respectively. This information will be written in a dedicated tracking 
table of the HFC. 

3.3. Filaments Behaviour 

In addition to the tracking algorithm, a module has been implemented in the 
code to characterize the behaviour of the filaments during their crossing the solar 
disk. To achieve this goal, the procedure checks, using the tracking identification 
numbers, that each co-rotating filament is well detected on every image between 
its first and last times of observation. If during this period the filament disappears 
(i.e., is not seen on one or more successive images), then the times corresponding 
to the last observation before its disappearance and the first observation after its 
re-appearance are saved. Moreover, if the filament appears once it has crossed 
the east limb, and/or disappears before reaching the west limb, the information 
is also returned by the code. To proceed in this case, average heliographic coordi- 
nates {<(ph>,<Xh> of the filament skeleton centres on the first/last image are 
used to calculate the predicted longitudes on the observation just before/after; 
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if these longitudes are less than 70° (in absolute value), then we assume that 
there is an appearance/disappearance after/before the limb. As for the tracking 
data, the results of the analysis are stored in the tracking table of the HFC. 

4. Performances of the Algorithm 
4.1. Assessments 

To optimize the tracking code, we need to define the combination of input 
parameters [aq, 9o,lo, Aio] that gives the best matching rate. This set will serve 
as a reference for tracking the filaments on the Meudon Ha spectroheliograms. 
At the same time, the assessment of these parameters will also give an indication 
of the performance of the algorithm, providing a positive tracking rate [Rp] after 
each run (as explained below). 

To find the best set of input parameters, we compared the code results with a 
representative sample of filaments manually tracked. This sample contains ap- 
proximately 2000 filaments, detected by the recognition code between 2000 and 
2009 on Meudon spectroheliograms. The tracking was then performed identifying 
"by-hand" the co-rotating filaments on synthetic synoptic maps, such as the one 
displayed on Figure[4j Direct checking on corresponding Ha images was also done 
to confirm the first choice. In order to keep a track of this selection, the results are 
written in ASCII format files using identification numbers (t'™''", T-man^ ^man-j 
in the automatic code. In practice, since all filaments were loaded from the HFC, 
the identification numbers [i/f^^^^ and Vi] are actually identical; this criterion will 
permit us to identify individual features between both sets of results. 
Once the sample was generated, we then proceeded as follows: 

i) The values of the input parameters ds, rmax: Wa, Wg, Wi, Pq, and were 
first fixed to a given value. (The effects of these parameters on the code 
performance were also tested, but the results of this evaluation are not pre- 
sented here to not overload the discussion.). We took respectively, ds = 1°, 
7'max = 5°, Wa = Wg = Wi = 1/3, Pq = Q = 0.5, and the ratio of minimum 
to maximum lengths of the compared filaments had to be greater than 1/3. 
Since the application of the tracking algorithm to too small filaments is neither 
relevant nor efficient, we thus decided to process only the filaments for which 
the length is equal or greater than 5°. 

ii) We carried out a series of executions of the code, using different combinations 
of values for the input parameters Oq, 9o, Iq, and Atg (reasonable ranges of 
values were previously defined for each parameter). After each computation, 
a positive tracking rate [Rp] was systematically calculated to evaluate the 
performance of the automated process compared to the manual one. In order 
to estimate Rp, we compare two by two the groups of filaments from both 
manual and automatic sets, which can be identified by their own tracking 
identification numbers r"^^" and respectively. For two filaments [i and j] 
belonging to each type of groups, if the condition i/,™^"=i/j is fulfilled, then 
Rp was incremented by 1/n where n is the total number of comparisons. 
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Table 1. Results of the assessments for the track- 
ing over one crossing on the solar disk: the positive 
rate Rp in percent (second line) and the corre- 
sponding number of filaments Nfn (third line), 
using the best input parameters, are indicated for 
different lengths of filament Zgi^c (first line). 





^skc 


> 5° 


> 6° 


> 8° 


> 10° 


> 15° 


Rp 
Ntii 


86% 
1991 


89% 
1698 


92% 
1263 


94% 
924 


95% 
479 



in) Finally, we keep the combination of values of aq, 9q, Iq, and A^o that gives 
the highest positive rate. 

In the first computation series, we did not take account of the results of the 
tracking over two rotations including Wj, then in a second series, we performed 
the same operation including both rf^^-n i^man compute Rp. Results of 
the computed series are presented in the following section. 

4.2. Results 

Table [l] shows the results of a first series of comparisons between the automated 
and the manual processes. As explained in the previous section, this first series 
does not take account of the tracking from one rotation to another using uf^^'^. 
The values of the input parameters given the highest rate Rp are in this case: 
ao = 5°, ^0 = 30°, and Iq =4°. Rp and the corresponding numbers of tracked 
filaments, which have a length of skeleton [Zgkc] greater than 5°, 6°, 8°, 10°, and 
15° respectively, are indicated. {Rp is provided in percent.) 
According to Table [T] the performance of the code increases when the length Zgkc 
increases: passing from 86% of successful tracking for the A^fii — 1991 filaments 
that have a length greater than 5°, to 95% for the Nfn — 479 filaments for which 
^ske ^ 15°- This is the consequence of the algorithm sensitivity to the length 
of the skeletons, since it is more difficult for the curve-matching to compare 
filaments that have too small lengths. 

Table [2] shows the results of the second series of assessment, including here 
the tracking results from a rotation to the following. In this case, the input 
parameters that give the best positive rate Rp are : ao = 5° , Oq — 20°, Iq = 2° , 
and A^o = 112 decimal days. The total positive rate [Rp] (given in percent) and 
the corresponding numbers of tracked filaments, which have a length of skeleton 
[/skc] greater than 5°, 6°, 8°, 10°, and 15° respectively, are provided. 
The method also offers good results for the tracking over successive rotations 
{i.e. 90% for ^skc ^ 15°), but the performance decreases more rapidly when /gkc 
decreases. We note that the best values of the input parameters are smaller 
than for the first series, because the algorithm tends to reduce the false tracking 
occurrences by using more restrictive values on larger features matched. 
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Table 2. Results of the assessments for two suc- 
cessive rotation tracking: the total positive rate 
Rp in percent (second line) and the correspond- 
ing number of filaments Nfn (third line), using 
the best input parameters for different lengths of 
filament Ic^y^^ (first line). 





^skc 


> 5° 


> 6° 


> 8° 


> 10° 


> 15° 


Rp 
Ntii 


74% 
1991 


77% 
1698 


81% 
1263 


88% 
924 


90% 
479 
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Figure 7. Same synoptic map as Figure |4] but after the execution of the tracking code using 
the best input parameters found in the first series. Each group of skeletons that actually 
correspond to the same co-rotating filament (i.e., which have the same tracking identification 
number) are plotted with one specific color. Only filaments that have a length greater than 
5° are tracked, and so plotted on the map. The two black dashed lines delimit the longitude 
range of the image taken at the time t, and for which tracking is initially performed. 



An example of a synoptic map, produced after a run of the tracking code 
using the best set of inputs parameters found in the first series, is shown on 
Figure [7] The colors of skeletons are proportional to the corresponding tracking 
identification numbers, in such a way that associated features appear with the 
same color on the plot. In most cases, we can notice that the larger features are 
well grouped by colors. 
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5. Conclusions 

We have presented here a new method to track filaments on solar images. This 
technique is based on a curve-matching algorithm, applied to the skeletons of 
filaments in a reference frame rotating with the Sun's surface. The comparison 
of the automated code results with a representative sample of filaments tracked 
manually confirms the good performance (close to w 90% ; depending on the 
length of the skeletons) of the method. Moreover, the results show the good 
aptitude of the process to identify the main parts of a segmented filament (as 
a complementary tool to the detection codes, that often offer such a process 
on each image). Application of the code to track features over two successive 
rotations gives also good results, but appears to be less efficient, especially for 
the smaller filaments; improvements to include smaller filaments cases are in 
progress. 

For now, a version of the code is already used to provide tracking information to 
the filaments table of the HFC. Stored content only concerns filaments detected 
on the Meudon Ha spectroheliograms, but efforts are currently made to also 
include BBSO observations. Joint use of these data sets will allow us to improve 
the ability of the code to follow filaments, by increasing the temporal resolution 
(Meudon Observatory provides only one image per day on average), but also 
will refine the analysis of the filament behaviour, since a low cadence limits the 
capability of the algorithm to accurately detect the occurrence time of events 
such as filament disappearances. 

In addition, several applications of this method are envisaged in the future. 
Firstly, automated tracking will be a primary step to the creation of filaments 



synthetic synoptic maps (Mouradian 1998 Aboudarham et al.\ 2007) ; these 



maps provides useful information about solar features to the community. Since 
the detection and tracking codes can be run independently (and insofar as a 
dedicated detection code can previously provide the pruned skeletons), inves- 
tigations are in progress to apply this technique to filament extraction codes 



working on other wavelengths observations (Buchlin et al. 2010). Finally, an 
extension to others solar features such as active regions or coronal holes, could 
be reasonably envisaged using appropriate matching algorithms. 
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